home *** CD-ROM | disk | FTP | other *** search
/ Nautilus 1992 July / Nautilus-3-8 / Nautilus-3-8.bin / Tools & Utilities / Techy Stuff / Source ƒ / sky source ƒ / MERC.C < prev    next >
Encoding:
C/C++ Source or Header  |  1992-06-16  |  3.2 KB  |  173 lines
  1. #include "sky.h"
  2.  
  3. extern struct merctab
  4. {
  5.     float f[2];
  6.     char c[3];
  7. } merctab[];
  8.  
  9. merc()
  10. {
  11.     double pturbl, pturbb, pturbr;
  12.     double lograd;
  13.     double dele, enom, vnom, nd, sl;
  14.     double q0, v0, t0, m0, j0 , s0;
  15.     double lsun, elong, ci, dlong;
  16.     double planp[7];
  17.     struct merctab *pp = &merctab[0];
  18.     double olong;
  19.     double temp;
  20.  
  21. /*
  22.  *    The arguments nnd coefficients are taken from
  23.  *    Simon Newcomb, Tables of the Heliocentric Motion
  24.  *    of Mercury
  25.  *    A.P.A.E. VI, part 2 (1895).
  26.  *
  27.  *    Here are the mean orbital elements.
  28.  */
  29.  
  30.     object = "Mercury     ";
  31.     ecc = .20561421 + .00002046*capt - 0.03e-6*capt2;
  32.     incl = 7.0028806 + .0018608*capt - 18.3e-6*capt2;
  33.     node = 47.145944 + 1.185208*capt + .0001739*capt2;
  34.     argp = 75.899697 + 1.555490*capt + .0002947*capt2;
  35.     mrad = .3870986;
  36.     anom = 102.279381 + 4.0923344364*eday + 6.7e-6*capt2;
  37.     motion = 4.0923770233;
  38.  
  39.     incl *= radian;
  40.     node *= radian;
  41.     argp *= radian;
  42.     anom = fmod(anom, 360.)*radian;
  43.     motion *= radian;
  44.  
  45. /*
  46.  *    Conventional mean anomalies of perturbing planets.
  47.  */
  48.  
  49.     q0 = 102.28  + 4.092334429*eday;
  50.     v0 = 212.536 + 1.602126105*eday;
  51.     t0 = -1.45  + .985604737*eday;
  52.     m0 = 319.66 + .524028480*eday;
  53.     j0 = 225.36 + .083086735*eday;
  54.     s0 = 175.68 + .033455441*eday;
  55.  
  56.     q0 *= radian;
  57.     v0 *= radian;
  58.     t0 *= radian;
  59.     m0 *= radian;
  60.     j0 *= radian;
  61.     s0 *= radian;
  62.  
  63.     planp[1] = q0;
  64.     planp[2] = v0;
  65.     planp[3] = t0;
  66.     planp[4] = m0;
  67.     planp[5] = j0;
  68.     planp[6] = s0;
  69.  
  70. /*
  71.  *    Computation of long period terms affecting the mean anomaly.
  72.  */
  73.  
  74.     anom = anom;
  75.  
  76. /*
  77.  *    Computation of elliptic orbit.
  78.  */
  79.  
  80.     enom = anom + ecc*sin(anom);
  81.     do {
  82.         dele = (anom - enom + ecc * sin(enom)) /
  83.             (1. - ecc*cos(enom));
  84.         enom += dele;
  85.     } while(fabs(dele) > 1.e-8);
  86.     vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
  87.             cos(enom/2.));
  88.     rad = mrad*(1. - ecc*cos(enom));
  89.  
  90. /*
  91.  *    Perturbations in longitude.
  92.  */
  93.  
  94.     pturbl = 0.;
  95.     for(;;){
  96.         if(pp->f[0]==0.){
  97.             pp++;
  98.             break;
  99.         }
  100.         pturbl += pp->f[0]*cos(pp->f[1] + pp->c[0]*q0 + pp->c[1]*planp[pp->c[2]]);
  101.         pp++;
  102.     }
  103.  
  104. /*
  105.  *    Perturbations in latitude.
  106.  */
  107.  
  108.     pturbb = 0.;
  109.     for(;;){
  110.         if(pp->f[0]==0.){
  111.             pp++;
  112.             break;
  113.         }
  114.         pturbb += pp->f[0]*cos(pp->f[1] + pp->c[0]*q0 + pp->c[1]*planp[pp->c[2]]);
  115.         pp++;
  116.     }
  117.  
  118. /*
  119.  *    Perturbations in log radius vector.
  120.  */
  121.  
  122.     pturbr = 0.;
  123.     for(;;){
  124.         if(pp->f[0]==0.){
  125.             pp++;
  126.             break;
  127.         }
  128.         pturbr += pp->f[0]*cos(pp->f[1] + pp->c[0]*q0 + pp->c[1]*planp[pp->c[2]]);
  129.         pp++;
  130.     }
  131.     pturbr *= 1.e-6;
  132.  
  133. /*
  134.  *    reduce to the ecliptic
  135.  */
  136.  
  137.     olong = vnom + argp + pturbl*radsec;
  138.     nd = olong - node;
  139.     lambda = node + atan2(sin(nd)*cos(incl), cos(nd));
  140.  
  141.     sl = sin(incl)*sin(nd);
  142.     beta = atan2(sl, sqrt(1.-sl*sl)) + pturbb*radsec;
  143.  
  144.     lograd = pturbr*2.30258509;
  145.     rad *= 1. + lograd;
  146.  
  147. /*
  148.  *    Compute motion for planetary aberration.
  149.  */
  150.  
  151.     temp = motion*mrad*mrad*sqrt(1.-ecc*ecc)/(rad*rad);
  152.     ldot = temp*sin(2.*(lambda-node))/sin(2.*(olong-node));
  153.     bdot = temp*sin(incl)*cos(lambda-node);
  154.     rdot = motion*mrad*ecc*sin(olong-argp)/sqrt(1.-ecc*ecc);
  155.  
  156. /*
  157.  *    Compute magnitude.
  158.  */
  159.  
  160.     lsun = 99.696678 + 0.9856473354*eday;
  161.     lsun *= radian;
  162.     elong = lambda - lsun;
  163.     ci = (rad - cos(elong))/sqrt(1. + rad*rad - 2.*rad*cos(elong));
  164.     dlong = atan2(sqrt(1.-ci*ci), ci)/radian;
  165.     mag = -.003 + .01815*dlong + .0001023*dlong*dlong;
  166.  
  167.     semi = 3.34;
  168.  
  169.     helio();
  170.     geo();
  171.  
  172. }
  173.